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Abstract 

We introduce a technique to solve numerically the relativistic Euler's equa- 
tions in scenarios with spherical symmetry using the standard Smoothed 
i ■ Particles Hydrodynamics method in cartesian coordinates. This implemen- 

(3jT)' tation allow us to increase the resolution of the simulations in order to obtain 

accurate results. We test our implementation studying the evolution of a per- 
fect fluid in a blast wave configuration in a fixed space-time . The technique 
can be easily generalized to axial symmetric problems. 
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1. Introduction 

Physical scenarios involving fluids are studied using different numerical 
methods. One of these standard methods is the Smoothed Particle Hydro- 
^ ; dynamics (SPH). 

Different implementations for Newtonian and relativistic Euler's equa- 
tions in three spatial dimension using SPH have been studied for many years 
0, 0, 0, 0, [fl] . The idea is to write the evolution equations in cartesian coordi- 
nates and using a Lagrangian scheme we follow the evolution of the elements 
of the fluid during the simulation. When the physical problem has spherical 
symmetry, the standard approach is to rewrite the evolution equations in 
spherical coordinates, try to find the best way to adjust the parameters of 



the discretization and then evolve the system under that symmetry [10 
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In this article we use the ideas introduced in |7| (in the context of nu- 
merical evolutions of black holes) in order to evolve the system of equations 
written in cartesian coordinates without rewriting the system of equations, 
only using the symmetries of the problem. This "Cartoon SPH" technique, 
can easily be generalized to systems with axial symmetry. It is very helpful 
and straightforward task if a standard SPH code is already working, pro- 
viding a simple way to obtain high resolution to evolve the systems with 
symmetries and obtain accurate results. 

The structure of the paper is the following: In section [2] we describe the 
standard SPH. In section [3] we describe relativistic Euler's equations. Then, 
in section H] we describe the idea and the implementation of the cartoon SPH. 
In section we present tests for the cartoon SPH implementation and finally, 
in section [6] we conclude. 



2. Standard SPH 

The Smoothed Particle Hydrodynamics is a method used to solve nu- 
merically hydrodynamical equations. It is a mesh free method, also called a 
Lagrangian method, because we are not dealing with a fixed grid, instead we 
use several nodes called particles distributed on the volume of the fluid that 
we are studying. 

The discretization of the functions and their derivatives in the SPH method, 



is carried out in two steps [14], |2| 



1. Integral representation of a function: Let / be a real valued function 
from R 3 — > R. We use the following identity: 

/(r) = / S(r - r')dr' « / f{r)W{r - r', h)dr' (1) 
J Ri Jn 

where the delta function has been approximated by the function W 
called the kernel and h -called the smoothing length of the kernel- 
defines the region where the kernel is different from zero, i.e. Q C R 3 . 
The kernel is a smooth function over R 3 specifically over Q and it is 
normalized to the unity according to J n Wdv = 1. W is assumed to be 
symmetric, i.e. it only depends on the norm of the vector r — r'. 

2. Particle approximation: We change the integration by a sum over dis- 
crete volume elements dr — > AV = l/n(r) where n is the number 
density, subdividing the fluid in N parts we get 
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N 



[ f(r)W(r - r', h)dr' « V ^W o6 :=< / > a , (2) 

where we use the convention that for any real valued function f a := 
f(r a ), with a = 1, . . . , N. 

The derivatives of any real valued function are obtained using the compact 
support 



V/(r) = / Vf(r)W(\\r-r'\\)dr> 
Jn 

= I f(r')VW(\\r-r'\\)dr' + I f(r)'W(\\r - r'\\)nds . (3) 
Jn J an 

The integral over the boundary of the volume dQ is zero because of the 
compact support of the kernel. Then 

< V / >«= J2 -VW ab . (4) 

Following the same procedure we obtain the approximation for the diver- 
gence of a vector: 

< dif > a = V - f b ■ V a W ab . (5) 

3. Hydrodynamic Relativistic Equation 

If we want to study the behavior of a fluid in a curved space time we need 
to use the laws of thermodynamics in curved-space-times, i.e., local baryon 
conservation, the first and second laws of thermodynamics plus the local law 
of energy-momentum conservation [161 ]: 

V • T = 0. (6) 

Choosing a coordinate basis {c^}, we express equation (jSj) as V M T MV = 0. 
Now, we assume the fluid can be approximated by a perfect fluid repre- 
sented by the following stress-energy tensor: T^ v = (pw + q)u ll u u + g >J,u (p + q). 
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Here p is the rest mass-energy density, w = l + e + p/pis the relativistic spe- 
cific enthalpy, e is specific internal energy, p is the hydro dynamic pressure 
and q is a quantity known as the artificial viscosity [H, 16 . 



In order to evolve the system for a Lagrangian formulation of relativistic 
hydrodynamic equations we need to do a 3 + 1 splitting of the space-time. 
The standard way to do this is using the ADM formalism, where the space- 
time is decomposed into an infinite foliation of spatial hyper-surfaces E t of 
constant t coordinate. The line element is given by 



ds 2 = g^dx^dx" = -(a 2 - f3 i f3 i )dt 2 + 2(3 i dx i dt + r) ij dx i dx j 



(7) 



where a is the lapse function, (3 l the shift vector, and 77^ the induced 3-metric 
on E t JjJ [l8|, |l7|. Greek indices run from to 3 and Latin indices from 1 to 
3. 

We can express the quantities either in the coordinate basis {d^} or in the 
basis formed by the Eulerian observer 4— velocity n and the spatial vector 
basis di, {n, di} (with n and di orthogonal). For example, we can express the 
4- velocity of a fluid in these two basis: u = 7(71 + v l di) or u = -{d t + v l di) 
where v % = av l — (3 l and 7 is the Lorentz factor. 

In order to obtain the Lagrangian equations of relativistic hydrodynamics 
we define the Lagrangian derivative as 



d ol 

— = -u^df, = d t + v l di , 
dt 7 

and we write the equations of motion in terms of this derivative. 
The local conservation of baryon number is given by 



(8) 



dD* 



+ D*d iV l = 0, 



(9) 



where 

D* = yf^g-ip/a = y/rj^p, (10) 

is the relativistic rest mass density, and y/fj = \J—g /a is the determinant of 
the induced 3— metric. 

With the spatial part (p = i) of equation (EJ) we obtain the relativistic 
momentum equation 



-S =- — 

dt 1 D* 



di [V=g(p + ?)] - ^T^d i9aP 



(11) 



4 



with the relativistic specific momentum Si defined by 



Si 



w H — 
P 



(12) 



Notice that equation f fTTj) contains spatial derivatives of the metric of the 
space-time. For simplicity, we assume that the fluid interacts very weakly 
with the space-time, and we are going to use a fixed background. This means 
that the metric is given during all our simulation and we can compute the 
spatial (and also the temporal) derivatives either analytically or numerically. 

The temporal part (// = 0) of equation give us the relativistic energy 
equation: 



dE 



1 

D~* 



with E = aE — (3 l St and E is the total relativistic specific energy 



E 



. 9 
w H — 

P 



1 ~ 



p + q 

PI 



(13) 



(14) 



Finally the system of equations must be closed with an equation of state 
P = p{p-> e )i we are g°i n g to use an ideal gas equation of state 



p = (r - l)pe. 



(15) 



3.1. Discretization of Motion Equations 

Using the ideas and equations introduced in section El the equation of 
the relativistic momentum can be written as: 



dt 



'9a ^2 mb 

b 



Pa + qab , Pb + qba 



+ 



V a W ab 



-9a 



D* 



(Pa + q a )V a (In v^) a - -TfV a [g aP ) a 



(16) 



where S a = {Si} a , and the metric gradients Vliiy/^ and Vg^ can be 
calculated from the given metric. 



5 



The relativistic energy equation is 



dE a 
dt 



E 



Pa + qab + Pb + qba 



D 



*2 



(v a + v b )-V a W ab 



D* 



(p a + q a )v a -V a (In v^) a + -Tf(g a ^ a 



• (17) 



Finally, there are two possible ways to obtain the density: The first one 
is recovering the density by summation, using equation (|2J) 



D * a = J2m b W ab 



(18) 



where we have used n a = D* a /m a . 

The second one is integrating the density using equation 

d 



—D* a = -^2m b (v b -v a )-V a W ab . 



(19) 



It can be proved that keeping h constant, these two equations are equiv- 
alent. 

3.2. Artificial Viscosity 

In order to handle the shocks that appear evolving Euler equations, we 
use the extra term q as we mentioned before. This term is the artificial 
viscosity and it is inspired in the standard artificial viscosity used in 0,0, 

In our simulations, we have used the following artificial viscous pressure 



q a = -Y, b m b 



q a b qba 



w, 



ab 



(20) 



q a b 



PaW a 




if (V -v)<0 



where 

-ac a h a (V -v) a + ph 2 a (V -v) 2 a 

otherwise 

The divergence of the velocity for the a— th particle is 

v a b-r a b 



(21) 



V • v) 



\r a b\ 2 + eh 2 ab 



(22) 
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where c a = \/Tp a / (p a w a ) is the relativistic sound velocity measured in the 
rest frame of the fluid, a, (3 and e are numerical parameters, v a b, r a b were 
define as the differences v a b = v a — v^, r a b = r a — r&, and h a b = {h a + hb)/2 
is the mean value of the smoothing lengths of particles a and b. This version 
of artificial viscosity [l| is equivalent to the invented by Monaghan et. al in 
0. 

It is also possible to use other methods that offer to solve numerically the 
discontinuities using a solution of the Riemman problem for the fluid motion 
equation [lil . fl3| . That kind of analysis will be part of future works. 



3.3. About the Implementation 

At the initial time to we have the initial data of our physical problem 
(r,v,p, p,e)t - With this information, we reconstruct the initial relativistic 
variables (D*, S,E) to using equations (1TU| [T2| [Tl"|) . We can now integrate 
the evolution equations to obtain the relativistic variables at the new time 
(D*,S,E) 

to+St- 

The next step is to recover the physical variables. This can be accom- 
plished solving numerically an algebraic equation for 7 

= (S 2 - E 2 ) 7 4 + 2GE^ + [E 2 - 2GS 2 - G 2 ) 7 2 - 2GE 7 + G 2 (l + S 2 ) 

(23) 

where 



(24) 



and 



E = E + — 
TD 



Once the value of 7 is known, it can be calculated the rest-mass density 
p using equation ( jTOj) . then the thermodynamic pressure p from the equation 

w + q / p =(E 1 -G)/( 1 2 -G) (25) 

(G = 1 — 1/r) and the equation of state ( |T5i) . also we can obtain the specific 
internal energy e from e = (f~jr; and finally the velocity from equation 
( fl~2l) using r^Sj = (u> + q/ p)^v l . 
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4. Cartoon SPH 



Now we are going to describe how the symmetry of the problem can be 
used to improve the numerical calculations. We are going to describe the 
method using spherical symmetry. The generalization to cylindrical symme- 
try is straightforward. 

Using the cartesian and the spherical coordinate vector basis k) and 
(f , 6, 0) respectively, we identify the axis k with the radial direction f. Then, 
we subdivide the sphere in N s shells. We impose that each shell has the same 
mass with the following relation: 



with Mt the total mass of the configuration and rf the inner boundary 
position of shell i (i — 1, . . . , N s ). To assign the position of the real particles 
or nodes we use 



where rf is the position of the i-th real particle. 

The third step consists in generate a set of new particles around each real 
particle. This can be accomplished building an sphere of radius hi ( smoothing 
length) and subdividing it in N v virtual particles. In order to use the 3D 
SPH code we need to assign values of the physical quantities for each virtual 
particle. By construction we know the position and volume of the virtual 
particles, then, the physical values required ( r v , v v ,p v , p v , e v ) can be assigned 
interpolating the values of the real particles. We can compute the mass of 
each one of the virtual particles using the simple relation m v = p v AV v . 
Finally, we can obtain the auxiliar variables (D* v , S v , E v ) and use the 3D 
SPH equations of motion to evolve the N s particles using for each one N v 
virtual particles. 

4-1- Constructing the Virtual Particles 

Given a physical particle located at rf 1 , we construct a sphere of radius 
hi containing N n physical particles inside of it. 

We use spherical coordinates and split the 6 and angles in N e and 
parts respectively. The radial coordinate is subdivided in N r + 1 parts. The 




(26) 




(27) 
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additional subdivision of the radial coordinate corresponds to the central 
physical particle. 

The positions of the boundaries of the volume elements with respect to 
the position of the physical particle, are obtained in the following way: 

• The radius hi is subdivided in N r + 1 shells, so the first division is 

r b = hi/(N r + l). (28) 

The rest are obtained demanding that all the shells have the same 
volume, using the recurrence relation 

r b a = K_ 1 + 7V- 1 ^_ 1 (l-(7V r + l)-3)] wi th a = l,...,N r (29) 

• 6 is subdivided in Ng equal parts 

e b „ = /3~, with p = l,...,N and a ^ 1. (30) 

Nq 

• is subdivided in equal parts 

0^ = 7-—, with 7 = l,...,iV^, and a^l. (31) 

We identify the position of the virtual particles with the geometrical cen- 
ters of these volume elements: 



r 

a 



~.b _i_ „h 
v _ la' a—1 



2 



a* - 2/3-1 71 (V) 
e ? ~ (32) 
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27-I 2tt 
2 'Nl 



The first radius is r\ = |Ar because the central particle has radius Ar 



2 

and its volume is 



^ V central = \* ^ (33) 



with Ar = 



The volume of each element constructed can be calculated by an analyt- 
ical expression 

2ir 

A^, 7 = - — (cos(6$) - cos(^_ 1 )) (r» - r^) . (34) 

4-1.1. Assigning values to the virtual particles 

Now we must assign physical values to each one of the virtual particles in 
order to evolve the system with the SPH algorithm in cartesian coordinates. 

The virtual sphere has been subdivided entirely in N v = (N r + 1) • N e ■ 
volume elements with a virtual particle at the center. Lets assign to each 
one of them a physical value related with the physical quantities. 

Instead of interpolating for the N v virtual particles, we use a subset of 
Nfan = (N r + 1) • No auxiliary particles lying on res-plane (cjf = 0). 



x 




Figure 1: Left. Location of the virtual particles in the rrz-plane (<j> v = 0). We call this set 
of points the fan. Given the radius r™^ we interpolate using the real particles. Right. We 
show how to assign values to the virtual particles copying the values of the particles in 
the fan, due to the one to one relation between the particles in the fan and the particles 
at each <f>". 

The coordinates assigned to each one of the particles in the fan are given 

by 
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fan v ■ nv 

V% = o 

zQ = rf + r-coB^. (35) 
The central particle coincides with the physical particle: 

P central Pi ' 

U central ~ U i) 
v,radial _ z 
V central ~ "i ' 

P central Pi ' 

(36) 



and we assign physical values to the virtual particles, interpolating the 
values of the real particles using the radial distance to the origin. 

Once the fan has been filled, we can copy for all the virtual particles in 
the sphere rotating around the z-axis, i.e., for all the see Figured! 

For the velocity, we interpolate the radial velocity w , then we recon- 
struct the cartesian components of the velocity such that we can introduce 
them in the SPH 3D algorithm. We use the unitary vector pointing from the 
origin to the particle 

-v v,radial ~ fm\ 

V a, Pn = e ",/3,7 ( 37 ) 

with | 
and 



r a,/3 



< Al = r^sin^cos^, 

yl Al = r^sm^sin^, (39) 
< A7 = rr + r-oos^. 

We identify each (a, /3, 7) with a virtual particle and the implementation 
in the 3D code is straightforward. 
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4-2. The Virtual Particle Approximation 

Before presenting the simulations obtained using the Cartoon SPH we 
verify that the construction of the virtual particles was made in a consistent 
way. We present four convergence tests and we present the results in Figure 

m 

• Volume approximation: 

We check that the volume of the virtual particles and the positions 
assigned to them are correct. We compare for a given sphere of radius 
h its analytical volume Vh = ^-h 3 and the approximated value of the 
volume using the volumes of the virtual particles (|34p 

< V >:= ZjAVj (40) 

where AVj is the volume of each one of the N v particles. Figure |2fa) 
shows the behavior of the relative error {E v = \ Vh ~< v> \ ) as function of 
the smoothing length. We can observe that the relative error is always 
close to the round-off error of the computer. 

• Normalization of the kernel: 

We verify that the relation J WdV = 1 is properly satisfied. For the 
i—th physical particle we have 

< W >:= ZjWijAVj « 1, (41) 

then, we compute the error E w = |1.0— < W > | increasing the 
number of virtual particles. We present the behavior in Figure Efb). 

• Derivative of the kernel: 

We verify that the approximation for the derivative J (r — r')VW(r — 
r')dr' = 1 is satisfied. The discretization of this equation is 

< DfW >:= S^r, - r ; ) V;ll '; ; Al } « 1 . (42) 

We compute the relative error E D . W := |1.0— < D^W > \ and obtain 
proper convergence for all the components of the derivative. We show 
the error of the derivative in the z direction in Figure (EJ^c)). 
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• Density 

The last convergence test compares the numerical approximation of a 
discontinuous profile of density with an initial profile (p exa ct) in a given 
region of the space: 

< p >i-.= ZjPjWijAVj , (43) 

where the subindex % labels the position r^. Again, the relative error is 
E i := |Pe3: '" :t, '~ <p> ' 1 , and the results are presented in the figure flSTd)). 

Pexact,i 

It is important to notice that if we choose a fixed number of particles and 
change the value of the smoothing length, the relative error remains constant. 

Now we can proceed to present the main simulations used to test the 
cartoon implementation. 



5. Blast wave 

We are going to consider a spherical fluid distribution with two regions: 
the first one n G [0, R in ] with pi = 1.0 x 10 5 , e/ = 2.5 x 10" 5 , Uj = and the 
second one r n G [R in ,R T } with p H = 0.125 x 10 5 ,e// = 2.0 x 10" 5 ,t>// = 0, 
here Ri n = 50 and Rt = 100. We are using units where c = 1. 

We assume a flat spacetime properly described by the Minkowsky metric 
9/j.u — diag(—l, 1, 1, 1). It is clear that the lapse function is a = 1 and the 
shift vector ft = 0, j = 1, 2, 3. The 3 — metric induced on the hyper surface 
of constant time t, S t ,is 77^ = Sij (the Kronecker-delta). 

The parameters for the artificial viscosity are a = 1.0 , (3 = 2.0 and 
e = 0.01. 

To find the initial distribution of the real particles we use equation ( I2UI) 

((rt) 3 ~ (rU 3 ) = ^f, (44) 

with M infiut the mass contained in the inner and outer regions. In region 
/ the index i takes the values i = 1, . . . , N in and in region //, i = N in + 
1, . . . , N s = N in + N out . 

We consider r$ = as the origin and r s Ns = Rt the radius of the complete 
sphere. Then we can get easily the recurrence equation 
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• Region I: 



r *=(t4 + (r - )31 ' (45) 



1/3 



where i = 1, . . . , N in . 



• Region II: 



XU 3 +(^i) 3 V /3 , (46) 



where i = N in + 1, . . . , N s . 
The implementation of equation ( 12"7|) is straightforward 



c ( Ki) 3 +(^) 3 y /3 (47) 



2 

In this test we have used N in = N out = 400 and N s = 800 



5.1. The evolution 

We present in the Figure ([3]) the physical quantities for the blast wave con- 
figuration, with the parameters mentioned above. The figure contains three 
different columns corresponding to three different times t = 1200, t = 3200 
and t = 5000. In this panels we can appreciate the behavior of the physical 
quantities where the velocities of the particles are very small compared with 
the speed of light (v « c) (lo| . 

We observe (from left to right in Figure [3]) the existence of the head, 
rarefaction wave, tail, contact discontinuity and shock waves. 

The behavior of the simulation is similar to the shock tube with some 
differences in the profile of the velocity between the regions of rarefaction 
and shock. The boundaries of the rarefaction zone are called the head and 
the tail, d3,[l|]. 

The pressure is continuos in the zone between the tail and shock point (in 
the classical problem this is supported by the Rankine-Hugoniot conditions). 

In the specific internal energy we notice a difference between the shock 
tube and the blast wave, in the region after the rarefaction and the contact 
discontinuity. All these deformations are result of the spherical symmetry of 
the problem. 
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It is clear that all the evolutions have oscillations in the area of the con- 
tact discontinuity. This is the result of the implementation of the artificial 
viscosity. It is important to mention that if we increase the number of real 
particles in the tests this oscillations decrease as we present in Figure (jlj) and 
it is because the linear interpolation is better with more subdivision in the 
radial direction. We used three different values of N s = 200,400,800. It is 
also possible to used better interpolations to decrease the oscillations. 



6. Conclusions 

The implementation introduced in this article is a successful method that 
can be used to deal with problems in spherical symmetry with 3 dimensional 
codes instead of rewriting the equations in that symmetry. This implemen- 
tation is simple and can be also applied to problems with axial symmetry. In 
the same way, it can also be implemented for the Newtonian Euler equations 



15j. The tests we presented here, give a clear idea of the behavior of the 
cartoon and can be used to perform further analysis and studies of physical 
and astrophysical scenarios. 
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Figure 2: Panel (a) shows the error Ey (defined in the main text) as a function of the 
smoothing length. Each line is obtained with a different numbers of virtual particles. The 
parameters used for this analysis are N v — N r ■ Ng ■ N$, where N r = 2" • 10, Ng = 80 and 
= 80, and n = 0, 1, . . . , 6. Panel (b) shows the error E\y for the kernel approximation 
using the same number of virtual particles (N v ) as before. It is clear the convergence of 
the error to zero. Panel (c) shows the error E^w of the derivative of the kernel in the 
z direction. We observe the same behavior as in panel (b). Panel (d) presents the error 
Ep from a discontinuous initial value of the density profile. Comparing the exact and 
the approximated profiles, we notice that near the discontinuity, the relative error does 
not converge to zero, but in all the other regions it has a clear convergence for different 
numbers of virtual particles. 
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Figure 3: In this figure we present three columns: the first one corresponds to time 
t = 1200, the second to time t = 3200 and the third to time t = 5000. The rows are such 
that the first one is the density profile, the second is the specific internal energy, the third 
is the thermodynamical pressure and the fourth one is the velocity in the z direction, which 
corresponds to the radial velocity. All the physical quantities present the same structure 
as in the shock tube (the order of the head, tail, contact discontinuity and shock point), 
but with differences in the profile. 
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Figure 4: Panel (a) In this panel we present the oscillations around the contact discontinu- 
ity and how they decrease when we increase the number of real particles N s = 200, 400, 800. 
Panel (b) we can see how the numerical solution approximate the step function better if 
the number of particles increase as above. 
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